rm(list = ls())
library(pacman)
p_load(lme4,lmerTest ,performance, tidyverse,lubridate, lattice, devtools, VGAM, ggpubr, plotly)
library("knitr")
name <- data0.anal %>% distinct(iso_code)
# JOR is inconsistent data
hist(data0.anal$week.case.eff)

check <- data0.anal %>% filter(iso_code == "JOR") %>% select(iso_code, time, iso_code, total.cases, week.case.eff)

summary(data0.anal$school.close.red)
##   0   2   3 
##  68  61 471
summary(data0.anal$work.close.red)
##   0   2   3 
## 156 313 131
summary(data0.anal$restric.gather.red)
##   0   1   3   4 
## 139  73  53 335
setwd("F:/STUDY HCA IN TAIWAN/0. ideas for thesis/COVID19 and NPIs/Data analysis/dataset/data_sep")

week20 <- data0.anal %>% filter(time == "19") %>% select(country, population, total.cases)
## Adding missing grouping variables: `iso_code`
write.table(week20, file = "tt_case_2020_w20.csv", sep = ",")
check.time <- data0.anal %>% select(id, week.case.eff, growth.rate.manual, timeStart, timeEff)
## Adding missing grouping variables: `iso_code`
 # Group plot
plot.group.case <- ggplot(data0.anal, aes(time+2, week.case.eff)) +
  stat_smooth(method = "loess", se = F, size = 2) +
  stat_summary(aes(col = country), fun.y = mean, geom = "line", alpha = 0.5) +
  xlab("Week") +
  ylab(" Average weekly growth rate")+ 
  theme_bw() + 
  theme(strip.background = element_blank()) +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +
  theme(legend.position = "bottom") + 
  scale_x_continuous(breaks = seq(2, 25, by = 1))
## Warning: `fun.y` is deprecated. Use `fun` instead.
plot.group.case
## `geom_smooth()` using formula 'y ~ x'

plot.group.string <- ggplot(data0.anal, aes(time, string.idx)) +
  stat_smooth(method = "loess", se = F, size = 2) +
  stat_summary(aes(col = country), fun.y = mean, geom = "line", alpha = 0.5) +
  xlab("Week") +
  ylab(" Stringency index")+ 
  theme_bw() + 
  theme(strip.background = element_blank()) +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) + 
  theme(legend.position = "bottom") +
  scale_x_continuous(breaks = seq(0, 23, by = 1))
## Warning: `fun.y` is deprecated. Use `fun` instead.
plot.group.string 
## `geom_smooth()` using formula 'y ~ x'

cairo_pdf(filename = "F:/STUDY HCA IN TAIWAN/0. ideas for thesis/COVID19 and NPIs/Data analysis/results/plot_sep/fig1_NPI2020_Sep29.pdf",
          width = 18, height = 8, #inch
          onefile = TRUE, family = "Segoe UI")

merge1 <- ggarrange(plot.group.string, plot.group.case, hjust = -1,
          ncol = 2, nrow = 1)
## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'
merge1

dev.off()
## png 
##   2
ggplotly(plot.group.case)
## `geom_smooth()` using formula 'y ~ x'
ggplotly(plot.group.string)
## `geom_smooth()` using formula 'y ~ x'
#https://www.guru99.com/r-scatter-plot-ggplot2.html#2
# clean daily grow rate = 0, and NAs
## NAs = 36 obs 
## zero growth rate = 58

## Check zero growth rate
data0.anal1 <- data0.anal
check <- data0.anal1 %>% filter(week.case.eff == 0 )  
check1 <- data0.anal1 %>% filter(week.case.eff != 0 ) %>% select(id, week.case.eff) 
## Adding missing grouping variables: `iso_code`
## Check NAs
summary(data0.anal$week.case.eff)
##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
## 0.000000 0.004217 0.017600 0.039388 0.051100 0.609670
#probit transformation
data0.anal <- data0.anal %>% filter(week.case.eff != 0) %>% mutate(probit.case = probitlink( week.case.eff))

qqnorm(data0.anal$probit.case)

hist(data0.anal$probit.case)

boxplot(data0.anal$probit.case)

summary(data0.anal$test.policy.red)
## Warning: Unknown or uninitialised column: `test.policy.red`.
## Length  Class   Mode 
##      0   NULL   NULL
data0.anal <- data0.anal %>% mutate(test.policy.red = case_when(
  test.policy==0~0,
  test.policy==1~0,
  test.policy==2~1, 
  test.policy==3~2
))
data0.anal$test.policy.red <- factor(data0.anal$test.policy.red)

summary(data0.anal$inter.travel)
##   0   1   2   3   4 
##  22   9  37 202 270
data0.anal <- data0.anal %>% mutate(inter.travel.red = case_when(
  inter.travel==0~0,
  inter.travel==1~0,
  inter.travel==2~0,
  inter.travel==3~1,
  inter.travel==4~2
))
data0.anal$inter.travel.red <- factor(data0.anal$inter.travel.red)


summary(data0.anal$wear.mask)
##   0   1   2   3   4 
## 188  58 101  90 103
data0.anal <- data0.anal %>% mutate(wear.mask.red = case_when(
  wear.mask==0~0,
  wear.mask==1~1,
  wear.mask==2~1,
  wear.mask==3~2,
  wear.mask==4~3
))
data0.anal$wear.mask.red <- factor(data0.anal$wear.mask.red)

data0.anal <- data0.anal %>% mutate(wear.mask.red1 = case_when(
  wear.mask==0~0,
  wear.mask==1~0,
  wear.mask==2~1,
  wear.mask==3~2,
  wear.mask==4~3
))
data0.anal$wear.mask.red1 <- factor(data0.anal$wear.mask.red1)


data0.anal$pop_log <- log(data0.anal$population)
hist(data0.anal$pop_log)

data0.anal$pop_denlog <- log(data0.anal$pop_density)
hist(data0.anal$pop_denlog)

data0.anal$median_age.log <- log(data0.anal$median_age)
hist(data0.anal$median_age.log)

data0.anal$gdp.log <- log(data0.anal$gdp)
hist(data0.anal$gdp.log)

data0.anal$incomegroup1 <- factor(data0.anal$incomegroup, levels = c("Low income", "Lower middle income", "Upper middle income" , "High income"))
data0.anal <- data0.anal %>% mutate(incomegroup2 = case_when(
  incomegroup1 == "Low income" ~ "LMIC",
  incomegroup1 == "Lower middle income" ~ "LMIC",
  incomegroup1 == "Upper middle income" ~ "upper_income",
  incomegroup1 == "High income" ~ "high_income"
)) 
data0.anal$incomegroup2 <- factor(data0.anal$incomegroup2, levels = c( "LMIC", "upper_income","high_income"))
library(performance)
aic <- function(model){
  fix2 <- performance(model)
  x<- as.character(substitute(model))
  fix2 <- fix2 %>% 
    mutate( npi = x)
  }
model.school.cls <- lmer(probit.case ~ school.close.red +
                           time + (time|iso_code),
                     data= data0.anal, REML = F, lmerControl(optimizer = "bobyqa") )
summary(model.school.cls)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: probit.case ~ school.close.red + time + (time | iso_code)
##    Data: data0.anal
## Control: lmerControl(optimizer = "bobyqa")
## 
##      AIC      BIC   logLik deviance df.resid 
##    484.9    519.3   -234.5    468.9      532 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -4.1656 -0.5395 -0.0463  0.4680  4.5459 
## 
## Random effects:
##  Groups   Name        Variance  Std.Dev. Corr 
##  iso_code (Intercept) 0.2044411 0.45215       
##           time        0.0006648 0.02578  -0.42
##  Residual             0.1066855 0.32663       
## Number of obs: 540, groups:  iso_code, 30
## 
## Fixed effects:
##                     Estimate Std. Error         df t value Pr(>|t|)    
## (Intercept)        -1.230536   0.106689  57.448375 -11.534  < 2e-16 ***
## school.close.red2   0.054009   0.094487 517.578250   0.572  0.56784    
## school.close.red3  -0.211454   0.069518 497.657722  -3.042  0.00248 ** 
## time               -0.063180   0.005503  30.598704 -11.480 1.28e-12 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) sch..2 sch..3
## schl.cls.r2 -0.344              
## schl.cls.r3 -0.558  0.662       
## time        -0.385 -0.171 -0.035
x1 <- aic(model.school.cls)


model.work.cls <- lmer(probit.case ~ work.close.red +
                           time + (time|iso_code),
                     data= data0.anal, REML = F, lmerControl(optimizer = "bobyqa") )
summary(model.work.cls)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: probit.case ~ work.close.red + time + (time | iso_code)
##    Data: data0.anal
## Control: lmerControl(optimizer = "bobyqa")
## 
##      AIC      BIC   logLik deviance df.resid 
##    488.9    523.3   -236.5    472.9      532 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -4.1447 -0.5262 -0.0373  0.4362  4.5904 
## 
## Random effects:
##  Groups   Name        Variance  Std.Dev. Corr 
##  iso_code (Intercept) 0.2276444 0.47712       
##           time        0.0007713 0.02777  -0.53
##  Residual             0.1070298 0.32715       
## Number of obs: 540, groups:  iso_code, 30
## 
## Fixed effects:
##                   Estimate Std. Error         df t value Pr(>|t|)    
## (Intercept)      -1.299756   0.098404  35.768193 -13.208 2.53e-15 ***
## work.close.red2  -0.156166   0.057166 399.519932  -2.732  0.00658 ** 
## work.close.red3  -0.215738   0.053792 509.475989  -4.011 6.96e-05 ***
## time             -0.059369   0.005857  30.094125 -10.136 3.23e-11 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) wrk..2 wrk..3
## wrk.cls.rd2 -0.287              
## wrk.cls.rd3 -0.310  0.656       
## time        -0.502 -0.173 -0.023
x2<- aic(model.work.cls)



model.pubevent <- lmer(probit.case ~ cancel.pubevent.red +
                           time + (time|iso_code),
                     data= data0.anal, REML = F, lmerControl(optimizer = "bobyqa") )
summary(model.pubevent)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: probit.case ~ cancel.pubevent.red + time + (time | iso_code)
##    Data: data0.anal
## Control: lmerControl(optimizer = "bobyqa")
## 
##      AIC      BIC   logLik deviance df.resid 
##    494.6    524.7   -240.3    480.6      533 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.9901 -0.5414 -0.0326  0.4599  4.3079 
## 
## Random effects:
##  Groups   Name        Variance  Std.Dev. Corr 
##  iso_code (Intercept) 0.2179870 0.46689       
##           time        0.0006557 0.02561  -0.48
##  Residual             0.1094033 0.33076       
## Number of obs: 540, groups:  iso_code, 30
## 
## Fixed effects:
##                        Estimate Std. Error         df t value Pr(>|t|)    
## (Intercept)           -1.302306   0.100308  40.951960 -12.983 4.19e-16 ***
## cancel.pubevent.red2  -0.163539   0.056440 477.439874  -2.898  0.00393 ** 
## time                  -0.059070   0.005393  28.857875 -10.953 8.56e-12 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) cnc..2
## cncl.pbvn.2 -0.417       
## time        -0.461 -0.064
x3<- aic(model.pubevent)


model.restr.gather <- lmer(probit.case ~ restric.gather.red +
                           time + (time|iso_code),
                     data= data0.anal, REML = F, lmerControl(optimizer = "bobyqa") )
summary(model.restr.gather)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: probit.case ~ restric.gather.red + time + (time | iso_code)
##    Data: data0.anal
## Control: lmerControl(optimizer = "bobyqa")
## 
##      AIC      BIC   logLik deviance df.resid 
##    492.6    531.3   -237.3    474.6      531 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.9120 -0.5585 -0.0522  0.4664  4.6081 
## 
## Random effects:
##  Groups   Name        Variance  Std.Dev. Corr 
##  iso_code (Intercept) 0.2244701 0.47378       
##           time        0.0006785 0.02605  -0.44
##  Residual             0.1074566 0.32781       
## Number of obs: 540, groups:  iso_code, 30
## 
## Fixed effects:
##                       Estimate Std. Error         df t value Pr(>|t|)    
## (Intercept)          -1.302949   0.098989  37.681114 -13.163 1.14e-15 ***
## restric.gather.red1  -0.035307   0.083952 524.739493  -0.421 0.674243    
## restric.gather.red3  -0.147442   0.088206 500.481696  -1.672 0.095234 .  
## restric.gather.red4  -0.222364   0.060163 522.223255  -3.696 0.000242 ***
## time                 -0.057899   0.005508  29.903788 -10.512 1.46e-11 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) rst..1 rst..3 rst..4
## rstrc.gth.1 -0.204                     
## rstrc.gth.3 -0.238  0.406              
## rstrc.gth.4 -0.353  0.384  0.554       
## time        -0.417 -0.108 -0.128 -0.128
x4 <- aic(model.restr.gather)



model.cls.trans <- lmer(probit.case ~ close.transport.red +
                           time + (time|iso_code),
                     data= data0.anal, REML = F, lmerControl(optimizer = "bobyqa") )
summary(model.cls.trans)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: probit.case ~ close.transport.red + time + (time | iso_code)
##    Data: data0.anal
## Control: lmerControl(optimizer = "bobyqa")
## 
##      AIC      BIC   logLik deviance df.resid 
##    499.0    529.0   -242.5    485.0      533 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.8582 -0.5388 -0.0419  0.4456  4.6640 
## 
## Random effects:
##  Groups   Name        Variance  Std.Dev. Corr 
##  iso_code (Intercept) 0.2474150 0.49741       
##           time        0.0007383 0.02717  -0.55
##  Residual             0.1096458 0.33113       
## Number of obs: 540, groups:  iso_code, 30
## 
## Fixed effects:
##                        Estimate Std. Error         df t value Pr(>|t|)    
## (Intercept)           -1.368495   0.100817  34.658429  -13.57 1.96e-15 ***
## close.transport.red1  -0.093965   0.049971 514.440037   -1.88   0.0606 .  
## time                  -0.059710   0.005648  28.885222  -10.57 1.93e-11 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) cls..1
## cls.trnsp.1 -0.288       
## time        -0.559 -0.038
x5 <- aic(model.cls.trans)


model.stay.home <- lmer(probit.case ~ stay.home.red +
                           time + (time|iso_code),
                     data= data0.anal, REML = F, lmerControl(optimizer = "bobyqa") )
summary(model.stay.home)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: probit.case ~ stay.home.red + time + (time | iso_code)
##    Data: data0.anal
## Control: lmerControl(optimizer = "bobyqa")
## 
##      AIC      BIC   logLik deviance df.resid 
##    492.5    522.5   -239.2    478.5      533 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -4.0139 -0.5302 -0.0375  0.4581  4.5679 
## 
## Random effects:
##  Groups   Name        Variance  Std.Dev. Corr 
##  iso_code (Intercept) 0.2459927 0.49598       
##           time        0.0007357 0.02712  -0.51
##  Residual             0.1078785 0.32845       
## Number of obs: 540, groups:  iso_code, 30
## 
## Fixed effects:
##                  Estimate Std. Error         df t value Pr(>|t|)    
## (Intercept)     -1.366991   0.097769  31.340528 -13.982 5.13e-15 ***
## stay.home.red1  -0.150050   0.046519 533.691032  -3.226  0.00133 ** 
## time            -0.059139   0.005633  29.064622 -10.498 2.11e-11 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) sty..1
## stay.hm.rd1 -0.178       
## time        -0.543 -0.053
x6 <- aic(model.stay.home)


model.restr.move <- lmer(probit.case ~ restric.move.red +
                           time + (time|iso_code),
                     data= data0.anal, REML = F, lmerControl(optimizer = "bobyqa") )
summary(model.restr.move)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: probit.case ~ restric.move.red + time + (time | iso_code)
##    Data: data0.anal
## Control: lmerControl(optimizer = "bobyqa")
## 
##      AIC      BIC   logLik deviance df.resid 
##    497.0    527.0   -241.5    483.0      533 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -4.0133 -0.5325 -0.0323  0.4478  4.5814 
## 
## Random effects:
##  Groups   Name        Variance  Std.Dev. Corr 
##  iso_code (Intercept) 0.2481191 0.49812       
##           time        0.0007701 0.02775  -0.57
##  Residual             0.1092041 0.33046       
## Number of obs: 540, groups:  iso_code, 30
## 
## Fixed effects:
##                     Estimate Std. Error         df t value Pr(>|t|)    
## (Intercept)        -1.371390   0.099162  32.374319  -13.83 3.88e-15 ***
## restric.move.red2  -0.103169   0.043910 532.547488   -2.35   0.0192 *  
## time               -0.059713   0.005741  28.756405  -10.40 2.95e-11 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) rst..2
## rstrc.mv.r2 -0.222       
## time        -0.592 -0.030
x7 <- aic(model.restr.move)


model.inter.trav <- lmer(probit.case ~ inter.travel.red +
                           time + (time|iso_code),
                     data= data0.anal, REML = F, lmerControl(optimizer = "bobyqa") )
summary(model.inter.trav)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: probit.case ~ inter.travel.red + time + (time | iso_code)
##    Data: data0.anal
## Control: lmerControl(optimizer = "bobyqa")
## 
##      AIC      BIC   logLik deviance df.resid 
##    482.4    516.7   -233.2    466.4      532 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.8349 -0.5442 -0.0098  0.4428  4.7752 
## 
## Random effects:
##  Groups   Name        Variance  Std.Dev. Corr 
##  iso_code (Intercept) 0.2593476 0.50926       
##           time        0.0006816 0.02611  -0.46
##  Residual             0.1048992 0.32388       
## Number of obs: 540, groups:  iso_code, 30
## 
## Fixed effects:
##                     Estimate Std. Error         df t value Pr(>|t|)    
## (Intercept)        -1.284909   0.112827  49.711369 -11.388 1.84e-15 ***
## inter.travel.red1  -0.050367   0.071978 497.830959  -0.700 0.484413    
## inter.travel.red2  -0.286159   0.075038 485.356616  -3.814 0.000155 ***
## time               -0.058036   0.005462  29.082018 -10.626 1.58e-11 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) int..1 int..2
## intr.trvl.1 -0.453              
## intr.trvl.2 -0.458  0.727       
## time        -0.419 -0.040 -0.080
x8 <- aic(model.inter.trav)

model.test <- lmer(probit.case ~ test.policy.red +
                           time + (time|iso_code),
                     data= data0.anal, REML = F, lmerControl(optimizer = "bobyqa") )
summary(model.test)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: probit.case ~ test.policy.red + time + (time | iso_code)
##    Data: data0.anal
## Control: lmerControl(optimizer = "bobyqa")
## 
##      AIC      BIC   logLik deviance df.resid 
##    495.2    529.5   -239.6    479.2      532 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.8244 -0.5376 -0.0278  0.4330  4.7334 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev. Corr 
##  iso_code (Intercept) 0.248623 0.49862       
##           time        0.000638 0.02526  -0.61
##  Residual             0.109687 0.33119       
## Number of obs: 540, groups:  iso_code, 30
## 
## Fixed effects:
##                   Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)       -1.38791    0.09984  32.96732 -13.901 2.42e-15 ***
## test.policy.red1  -0.08567    0.06877 411.79676  -1.246  0.21355    
## test.policy.red2  -0.25225    0.08171 385.83663  -3.087  0.00217 ** 
## time              -0.05424    0.00572  36.55832  -9.481 2.17e-11 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) tst..1 tst..2
## tst.plcy.r1 -0.247              
## tst.plcy.r2 -0.132  0.486       
## time        -0.511 -0.261 -0.345
x9 <- aic(model.test)

model.tracing <- lmer(probit.case ~ contact.tracing +
                           time + (time|iso_code),
                     data= data0.anal, REML = F, lmerControl(optimizer = "bobyqa") )
summary(model.tracing)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: probit.case ~ contact.tracing + time + (time | iso_code)
##    Data: data0.anal
## Control: lmerControl(optimizer = "bobyqa")
## 
##      AIC      BIC   logLik deviance df.resid 
##    493.2    527.6   -238.6    477.2      532 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.8719 -0.5219 -0.0381  0.4368  4.7427 
## 
## Random effects:
##  Groups   Name        Variance  Std.Dev. Corr 
##  iso_code (Intercept) 0.3193535 0.56511       
##           time        0.0007502 0.02739  -0.70
##  Residual             0.1080664 0.32873       
## Number of obs: 540, groups:  iso_code, 30
## 
## Fixed effects:
##                    Estimate Std. Error         df t value Pr(>|t|)    
## (Intercept)       -1.806154   0.252055 128.554319  -7.166 5.37e-11 ***
## contact.tracing1   0.591454   0.251570 210.954308   2.351   0.0196 *  
## contact.tracing2   0.293960   0.236264 187.365946   1.244   0.2150    
## time              -0.060189   0.005721  29.197895 -10.521 1.90e-11 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) cntc.1 cntc.2
## cntct.trcn1 -0.880              
## cntct.trcn2 -0.894  0.934       
## time        -0.211 -0.087 -0.109
x10 <- aic(model.tracing)



model.mask <- lmer(probit.case ~ wear.mask.red +
                           time + (time|iso_code),
                     data= data0.anal, REML = F, lmerControl(optimizer = "bobyqa") )
summary(model.mask)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: probit.case ~ wear.mask.red + time + (time | iso_code)
##    Data: data0.anal
## Control: lmerControl(optimizer = "bobyqa")
## 
##      AIC      BIC   logLik deviance df.resid 
##    482.8    521.4   -232.4    464.8      531 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -4.3769 -0.5310 -0.0397  0.4861  4.4528 
## 
## Random effects:
##  Groups   Name        Variance  Std.Dev. Corr 
##  iso_code (Intercept) 0.2190665 0.46805       
##           time        0.0006316 0.02513  -0.65
##  Residual             0.1078840 0.32846       
## Number of obs: 540, groups:  iso_code, 30
## 
## Fixed effects:
##                  Estimate Std. Error         df t value Pr(>|t|)    
## (Intercept)     -1.356786   0.092316  29.572995 -14.697 3.94e-15 ***
## wear.mask.red1  -0.343768   0.068662 493.118928  -5.007 7.72e-07 ***
## wear.mask.red2  -0.223732   0.085326 332.601163  -2.622  0.00914 ** 
## wear.mask.red3  -0.125477   0.071999 404.955827  -1.743  0.08214 .  
## time            -0.050290   0.005932  42.460894  -8.478 1.11e-10 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) wr.m.1 wr.m.2 wr.m.3
## wer.msk.rd1 -0.143                     
## wer.msk.rd2 -0.065  0.446              
## wer.msk.rd3 -0.042  0.285  0.300       
## time        -0.558 -0.299 -0.354 -0.344
x11 <- aic(model.mask)

model.info <- lmer(probit.case ~ info.campaign +
                           time + (time|iso_code),
                     data= data0.anal, REML = F, lmerControl(optimizer = "bobyqa") )
summary(model.info)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: probit.case ~ info.campaign + time + (time | iso_code)
##    Data: data0.anal
## Control: lmerControl(optimizer = "bobyqa")
## 
##      AIC      BIC   logLik deviance df.resid 
##    500.9    535.2   -242.5    484.9      532 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.8260 -0.5417 -0.0377  0.4489  4.6673 
## 
## Random effects:
##  Groups   Name        Variance  Std.Dev. Corr 
##  iso_code (Intercept) 0.2448265 0.49480       
##           time        0.0006834 0.02614  -0.56
##  Residual             0.1101943 0.33196       
## Number of obs: 540, groups:  iso_code, 30
## 
## Fixed effects:
##                  Estimate Std. Error         df t value Pr(>|t|)    
## (Intercept)     -1.609405   0.217299 280.481998  -7.406 1.53e-12 ***
## info.campaign1   0.428322   0.251262 495.424675   1.705   0.0889 .  
## info.campaign2   0.172464   0.201308 462.073047   0.857   0.3920    
## time            -0.059719   0.005522  29.434366 -10.815 9.03e-12 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) inf.c1 inf.c2
## info.cmpgn1 -0.745              
## info.cmpgn2 -0.895  0.795       
## time        -0.191 -0.017 -0.088
x12 <- aic(model.info)
#population: no
model.pop <- lmer(probit.case ~ pop_log +
                           time + (time|iso_code),
                     data= data0.anal, REML = F, lmerControl(optimizer = "bobyqa") )
summary(model.pop)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: probit.case ~ pop_log + time + (time | iso_code)
##    Data: data0.anal
## Control: lmerControl(optimizer = "bobyqa")
## 
##      AIC      BIC   logLik deviance df.resid 
##    501.3    531.4   -243.7    487.3      533 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.8339 -0.5429 -0.0345  0.4384  4.6547 
## 
## Random effects:
##  Groups   Name        Variance  Std.Dev. Corr 
##  iso_code (Intercept) 0.2516624 0.50166       
##           time        0.0007352 0.02711  -0.60
##  Residual             0.1106882 0.33270       
## Number of obs: 540, groups:  iso_code, 30
## 
## Fixed effects:
##              Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept) -0.686257   0.698982 31.227979  -0.982    0.334    
## pop_log     -0.042756   0.040183 30.539823  -1.064    0.296    
## time        -0.060179   0.005643 28.806722 -10.664 1.63e-11 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##         (Intr) pop_lg
## pop_log -0.990       
## time    -0.097  0.008
y1 <- aic(model.pop)

#population density

model.popden <- lmer(probit.case ~ pop_denlog +
                           time + (time|iso_code),
                     data= data0.anal, REML = F, lmerControl(optimizer = "bobyqa") )
summary(model.popden)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: probit.case ~ pop_denlog + time + (time | iso_code)
##    Data: data0.anal
## Control: lmerControl(optimizer = "bobyqa")
## 
##      AIC      BIC   logLik deviance df.resid 
##    501.5    531.6   -243.8    487.5      533 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.8336 -0.5381 -0.0393  0.4384  4.6406 
## 
## Random effects:
##  Groups   Name        Variance  Std.Dev. Corr 
##  iso_code (Intercept) 0.2391485 0.48903       
##           time        0.0007352 0.02711  -0.57
##  Residual             0.1106907 0.33270       
## Number of obs: 540, groups:  iso_code, 30
## 
## Fixed effects:
##              Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept) -1.652227   0.252408 33.689187  -6.546 1.77e-07 ***
## pop_denlog   0.045603   0.046536 30.189707   0.980    0.335    
## time        -0.060053   0.005641 28.813550 -10.646 1.69e-11 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##            (Intr) pp_dnl
## pop_denlog -0.926       
## time       -0.240  0.009
y2 <- aic(model.popden)

# median age
model.age <- lmer(probit.case ~ median_age +
                           time + (time|iso_code),
                     data= data0.anal, REML = F, lmerControl(optimizer = "bobyqa") )
summary(model.age)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: probit.case ~ median_age + time + (time | iso_code)
##    Data: data0.anal
## Control: lmerControl(optimizer = "bobyqa")
## 
##      AIC      BIC   logLik deviance df.resid 
##    495.7    525.8   -240.9    481.7      533 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.8248 -0.5395 -0.0350  0.4409  4.6803 
## 
## Random effects:
##  Groups   Name        Variance  Std.Dev. Corr 
##  iso_code (Intercept) 0.2602538 0.5102        
##           time        0.0007343 0.0271   -0.71
##  Residual             0.1107129 0.3327        
## Number of obs: 540, groups:  iso_code, 30
## 
## Fixed effects:
##              Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept) -0.477093   0.341241 34.144716  -1.398  0.17110    
## median_age  -0.030262   0.010459 30.882357  -2.893  0.00693 ** 
## time        -0.060290   0.005645 28.726092 -10.681 1.62e-11 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##            (Intr) medn_g
## median_age -0.957       
## time       -0.217  0.007
y3 <- aic(model.age)


# gdp
model.gdp <- lmer(probit.case ~ gdp.log +
                           time + (time|iso_code),
                     data= data0.anal, REML = F, lmerControl(optimizer = "bobyqa") )
summary(model.gdp)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: probit.case ~ gdp.log + time + (time | iso_code)
##    Data: data0.anal
## Control: lmerControl(optimizer = "bobyqa")
## 
##      AIC      BIC   logLik deviance df.resid 
##      502      532     -244      488      533 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.8315 -0.5445 -0.0348  0.4372  4.6586 
## 
## Random effects:
##  Groups   Name        Variance  Std.Dev. Corr 
##  iso_code (Intercept) 0.2587356 0.5087        
##           time        0.0007347 0.0271   -0.61
##  Residual             0.1106916 0.3327        
## Number of obs: 540, groups:  iso_code, 30
## 
## Fixed effects:
##              Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept) -0.910145   0.733937 31.056382  -1.240    0.224    
## gdp.log     -0.053077   0.075327 30.370869  -0.705    0.486    
## time        -0.060187   0.005642 28.790591 -10.667 1.63e-11 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##         (Intr) gdp.lg
## gdp.log -0.991       
## time    -0.098  0.012
y4 <- aic(model.gdp)

# test per 1000 => too many missing because the early stage of pandemic
model.test1000 <- lmer(probit.case ~ test.thousand +
                           time + (time|iso_code),
                     data= data0.anal, REML = F, lmerControl(optimizer = "bobyqa") )
summary(model.test1000)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: probit.case ~ test.thousand + time + (time | iso_code)
##    Data: data0.anal
## Control: lmerControl(optimizer = "bobyqa")
## 
##      AIC      BIC   logLik deviance df.resid 
##    301.3    329.6   -143.7    287.3      411 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -4.3822 -0.4110 -0.0188  0.3532  5.1613 
## 
## Random effects:
##  Groups   Name        Variance  Std.Dev. Corr 
##  iso_code (Intercept) 0.1956395 0.44231       
##           time        0.0009157 0.03026  -0.64
##  Residual             0.0879862 0.29662       
## Number of obs: 418, groups:  iso_code, 24
## 
## Fixed effects:
##                 Estimate Std. Error         df t value Pr(>|t|)    
## (Intercept)   -1.3520967  0.0986303 23.2956823 -13.709 1.21e-12 ***
## test.thousand -0.0009576  0.0008348 44.7984633  -1.147    0.257    
## time          -0.0577606  0.0075758 26.9090064  -7.624 3.43e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) tst.th
## test.thosnd  0.099       
## time        -0.665 -0.414
y5 <- aic(model.test1000)


# percentage of urban
model.per_urban <- lmer(probit.case ~ per_urban + 
                          time + (time|iso_code), 
                        data= data0.anal, REML = F, lmerControl(optimizer = "bobyqa") )

summary(model.per_urban)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: probit.case ~ per_urban + time + (time | iso_code)
##    Data: data0.anal
## Control: lmerControl(optimizer = "bobyqa")
## 
##      AIC      BIC   logLik deviance df.resid 
##    502.0    532.1   -244.0    488.0      533 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.8322 -0.5443 -0.0350  0.4377  4.6549 
## 
## Random effects:
##  Groups   Name        Variance  Std.Dev. Corr 
##  iso_code (Intercept) 0.2530655 0.50306       
##           time        0.0007347 0.02711  -0.59
##  Residual             0.1106905 0.33270       
## Number of obs: 540, groups:  iso_code, 30
## 
## Fixed effects:
##              Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept) -1.307725   0.204141 36.795138  -6.406 1.82e-07 ***
## per_urban   -0.001889   0.002944 30.543970  -0.642    0.526    
## time        -0.060181   0.005642 28.816358 -10.667 1.61e-11 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##           (Intr) pr_rbn
## per_urban -0.878       
## time      -0.315  0.015
y6 <- aic(model.per_urban)

# uhc_index
model.uhc <- lmer(probit.case ~ uhc_index + 
                          time + (time|iso_code), 
                        data= data0.anal, REML = F, lmerControl(optimizer = "bobyqa") )

summary(model.uhc)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: probit.case ~ uhc_index + time + (time | iso_code)
##    Data: data0.anal
## Control: lmerControl(optimizer = "bobyqa")
## 
##      AIC      BIC   logLik deviance df.resid 
##    494.8    524.9   -240.4    480.8      533 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.8284 -0.5384 -0.0308  0.4440  4.6733 
## 
## Random effects:
##  Groups   Name        Variance  Std.Dev. Corr 
##  iso_code (Intercept) 0.2347485 0.48451       
##           time        0.0007333 0.02708  -0.68
##  Residual             0.1107152 0.33274       
## Number of obs: 540, groups:  iso_code, 30
## 
## Fixed effects:
##              Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept) -0.244745   0.404956 33.184394  -0.604  0.54970    
## uhc_index   -0.017419   0.005823 30.926225  -2.992  0.00541 ** 
## time        -0.060373   0.005639 28.815358 -10.706 1.48e-11 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##           (Intr) uhc_nd
## uhc_index -0.972       
## time      -0.178  0.015
y7 <- aic(model.uhc)

# health worker density
model.hworker <- lmer(probit.case ~ healthworker_den + 
                          time + (time|iso_code), 
                        data= data0.anal, REML = F, lmerControl(optimizer = "bobyqa") )

summary(model.hworker)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: probit.case ~ healthworker_den + time + (time | iso_code)
##    Data: data0.anal
## Control: lmerControl(optimizer = "bobyqa")
## 
##      AIC      BIC   logLik deviance df.resid 
##    499.2    529.3   -242.6    485.2      533 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.8323 -0.5420 -0.0319  0.4387  4.6618 
## 
## Random effects:
##  Groups   Name        Variance  Std.Dev. Corr 
##  iso_code (Intercept) 0.2443999 0.49437       
##           time        0.0007337 0.02709  -0.63
##  Residual             0.1107056 0.33272       
## Number of obs: 540, groups:  iso_code, 30
## 
## Fixed effects:
##                   Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)      -1.166817   0.168425 40.689289  -6.928 2.17e-08 ***
## healthworker_den -0.004002   0.002162 31.275638  -1.851   0.0737 .  
## time             -0.060325   0.005640 28.835170 -10.697 1.50e-11 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) hlthw_
## hlthwrkr_dn -0.821       
## time        -0.392  0.019
y8 <- aic(model.hworker)

#Income group
model.income <- lmer(probit.case ~ incomegroup2 + 
                          time + (time|iso_code), 
                        data= data0.anal, REML = F, lmerControl(optimizer = "bobyqa") )

summary(model.income)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: probit.case ~ incomegroup2 + time + (time | iso_code)
##    Data: data0.anal
## Control: lmerControl(optimizer = "bobyqa")
## 
##      AIC      BIC   logLik deviance df.resid 
##    500.2    534.6   -242.1    484.2      532 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.8123 -0.5392 -0.0420  0.4475  4.6250 
## 
## Random effects:
##  Groups   Name        Variance  Std.Dev. Corr 
##  iso_code (Intercept) 0.2109792 0.45932       
##           time        0.0007316 0.02705  -0.57
##  Residual             0.1107088 0.33273       
## Number of obs: 540, groups:  iso_code, 30
## 
## Fixed effects:
##                           Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)              -1.297804   0.123201 38.978726 -10.534 5.78e-13 ***
## incomegroup2upper_income -0.352801   0.168818 30.379305  -2.090   0.0451 *  
## incomegroup2high_income  -0.061070   0.174839 30.322594  -0.349   0.7293    
## time                     -0.060232   0.005628 28.965326 -10.703 1.40e-11 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) incmgrp2p_ incmgrp2h_
## incmgrp2pp_ -0.582                      
## incmgrp2hg_ -0.562  0.405               
## time        -0.462  0.018      0.017
y9 <- aic(model.income)

# SDI 
model.sdi <- lmer(probit.case ~ sdi_index + 
                          time + (time|iso_code), 
                        data= data0.anal, REML = F, lmerControl(optimizer = "bobyqa") )

summary(model.sdi)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: probit.case ~ sdi_index + time + (time | iso_code)
##    Data: data0.anal
## Control: lmerControl(optimizer = "bobyqa")
## 
##      AIC      BIC   logLik deviance df.resid 
##    501.5    531.6   -243.8    487.5      533 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.8302 -0.5449 -0.0345  0.4348  4.6649 
## 
## Random effects:
##  Groups   Name        Variance  Std.Dev. Corr 
##  iso_code (Intercept) 0.2622524 0.5121        
##           time        0.0007344 0.0271   -0.62
##  Residual             0.1106941 0.3327        
## Number of obs: 540, groups:  iso_code, 30
## 
## Fixed effects:
##              Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept) -1.060386   0.371535 33.263400  -2.854  0.00737 ** 
## sdi_index   -0.542060   0.536009 30.775256  -1.011  0.31976    
## time        -0.060230   0.005642 28.779814 -10.674 1.61e-11 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##           (Intr) sd_ndx
## sdi_index -0.964       
## time      -0.188  0.014
y10 <- aic(model.sdi)

# mobi.composite 
model.mobi <- lmer(probit.case ~ mobi.composite + 
                          time + (time|iso_code), 
                        data= data0.anal, REML = F, lmerControl(optimizer = "bobyqa") )

summary(model.mobi)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: probit.case ~ mobi.composite + time + (time | iso_code)
##    Data: data0.anal
## Control: lmerControl(optimizer = "bobyqa")
## 
##      AIC      BIC   logLik deviance df.resid 
##    393.6    422.9   -189.8    379.6      478 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -4.2626 -0.5368 -0.0019  0.4520  4.7071 
## 
## Random effects:
##  Groups   Name        Variance  Std.Dev. Corr 
##  iso_code (Intercept) 0.2106074 0.45892       
##           time        0.0005974 0.02444  -0.29
##  Residual             0.0977931 0.31272       
## Number of obs: 485, groups:  iso_code, 25
## 
## Fixed effects:
##                  Estimate Std. Error         df t value Pr(>|t|)    
## (Intercept)    -1.235e+00  1.024e-01  3.117e+01 -12.064 2.81e-13 ***
## mobi.composite  3.613e-03  8.774e-04  4.519e+02   4.118 4.54e-05 ***
## time           -6.415e-02  5.528e-03  2.454e+01 -11.604 1.88e-11 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) mb.cmp
## mobi.compst  0.341       
## time        -0.369 -0.072
y11 <- aic(model.mobi)

aic_table <- rbind(x1, x2, x3, x4, x6, x7, x8, x9, x10, x11, y3, y7)
model.mul <- lmer(probit.case ~ uhc_index  + inter.travel.red  + wear.mask.red +  test.policy.red  +

                                               time  + (time|iso_code),
                     data= data0.anal, REML = F, lmerControl(optimizer = "bobyqa") )
                                            summary(model.mul)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: probit.case ~ uhc_index + inter.travel.red + wear.mask.red +  
##     test.policy.red + time + (time | iso_code)
##    Data: data0.anal
## Control: lmerControl(optimizer = "bobyqa")
## 
##      AIC      BIC   logLik deviance df.resid 
##    464.5    524.6   -218.3    436.5      526 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -4.3417 -0.5247 -0.0151  0.4634  4.6497 
## 
## Random effects:
##  Groups   Name        Variance  Std.Dev. Corr 
##  iso_code (Intercept) 0.2186297 0.46758       
##           time        0.0005131 0.02265  -0.64
##  Residual             0.1025931 0.32030       
## Number of obs: 540, groups:  iso_code, 30
## 
## Fixed effects:
##                     Estimate Std. Error         df t value Pr(>|t|)    
## (Intercept)        -0.319755   0.423282  33.223032  -0.755  0.45532    
## uhc_index          -0.013109   0.006120  31.463177  -2.142  0.04003 *  
## inter.travel.red1  -0.059617   0.071211 494.020215  -0.837  0.40289    
## inter.travel.red2  -0.238736   0.074238 461.469385  -3.216  0.00139 ** 
## wear.mask.red1     -0.311671   0.068913 506.738186  -4.523 7.61e-06 ***
## wear.mask.red2     -0.184827   0.084831 333.205558  -2.179  0.03005 *  
## wear.mask.red3     -0.089169   0.070741 335.209161  -1.260  0.20837    
## test.policy.red1   -0.092283   0.067384 374.048339  -1.370  0.17166    
## test.policy.red2   -0.169078   0.081638 337.485494  -2.071  0.03911 *  
## time               -0.046106   0.005828  52.223350  -7.912 1.71e-10 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) uhc_nd int..1 int..2 wr.m.1 wr.m.2 wr.m.3 tst..1 tst..2
## uhc_index   -0.966                                                        
## intr.trvl.1 -0.068 -0.055                                                 
## intr.trvl.2 -0.128  0.013  0.732                                          
## wer.msk.rd1  0.073 -0.121  0.084 -0.003                                   
## wer.msk.rd2  0.076 -0.113  0.159  0.126  0.470                            
## wer.msk.rd3 -0.068  0.059  0.053  0.005  0.288  0.313                     
## tst.plcy.r1  0.008 -0.064  0.001 -0.100  0.163  0.017 -0.040              
## tst.plcy.r2  0.162 -0.175 -0.125 -0.196  0.000 -0.122 -0.151  0.502       
## time        -0.175  0.091 -0.063 -0.036 -0.327 -0.341 -0.299 -0.266 -0.263
performance(model.mul)
## # Indices of model performance
## 
## AIC     |     BIC | R2 (cond.) | R2 (marg.) |   ICC |  RMSE | Sigma
## -------------------------------------------------------------------
## 464.528 | 524.610 |      0.761 |      0.405 | 0.598 | 0.305 | 0.320
# model.mul1 <- lmer(probit.case ~
#  work.close.red + uhc_index + wear.mask +  test.policy  + incomegroup2 +
#                                                time  + (time|iso_code),
#                      data= data0.anal, REML = F, lmerControl(optimizer = "bobyqa") )
#                                             summary(model.mul1)
# performance(model.mul1)
# anova(model.mul, model.mul1)
p_load(broom.mixed)
sum_table <- tidy(model.mul, conf.int = T)

# create new table summary of model
new_table <- sum_table[2:15,]
 
term1 <- c(
"UHC Index\n",
"Stay at home\nwith minimal exception\n",
"Screening arrivals\n",
"Quarantine arrivals\nfrom some or all regions\n",
"Ban arrivals from\nsome regions\n",
"Ban on all regions\n",
"School close some \nlevels\n",
"School close all levels\n",
"Only those who both have\nsymptoms+ meet specific criteria\n",
"Testing of anyone \nshowing symptoms\n",
"Open public testing\n",
"Upper income countries\n",
"High income countries\n",
"Time\n"
  )
new_table1 <- cbind(new_table, term1)
new_table1$term2 <- factor(new_table1$term1, levels = new_table1$term1)

cairo_pdf(filename = "F:/STUDY HCA IN TAIWAN/0. ideas for thesis/COVID19 and NPIs/Data analysis/results/plot_sep/fig2_year2020_eff_NPI_Sep25.pdf",
          width = 15, height = 10, #inch
          onefile = TRUE, family = "Segoe UI")

p.year2020<- ggplot(new_table1, aes(x = term2, y = estimate,
                     ymin = conf.low, ymax = conf.high)) +
    geom_hline( yintercept = 0, color = 'Blue' ) + ylim(-1.5,1) + 
    geom_linerange() + geom_point() + coord_flip() + 
  ylab("Probit Transformation of Average Weekly Growth Rate ") + 
  xlab("Factors") +
  theme_minimal() + 
  theme(strip.background = element_rect(NA), 
  axis.text = element_text(size = 20),
  axis.title = element_text(size = 20),
  legend.text = element_text(size = 10),
  legend.title = element_text(size = 10),
  strip.text = element_text(size = 8),
  legend.position = "top")
p.year2020
## Warning: Removed 5 rows containing missing values (geom_segment).
## Warning: Removed 1 rows containing missing values (geom_point).
dev.off()
## png 
##   2
#devtools::install_github("goodekat/redres")
library(redres)

data0.anal1 <- data0.anal %>%  filter( is.na(uhc_index) ==F)

rc_resids <- compute_redres(model.mul)
pm_resids <- compute_redres(model.mul, type = "pearson_mar")
sc_resids <- compute_redres(model.mul, type = "std_cond")
resids <- data.frame(data0.anal1$iso_code, rc_resids, pm_resids, sc_resids)
plot_redres(model.mul, type = "std_cond")

plot_resqq(model.mul)

plot_ranef(model.mul)

plot(data0.anal1$time, rc_resids,ylim=c(-3,3))
abline(h=0, col="blue")

plot(data0.anal1$iso_code, sc_resids)
abline(h=c(0,2.96,-2.96), col="blue")

# # run other diagnosis
# p_load(patchwork, randomForest)
# 
# check_model(model.mul) 
# 
 plot(check_distribution(model.mul))
## Warning: Removed 4 rows containing missing values (geom_segment).
## Warning: Removed 4 rows containing missing values (geom_point).

# 
 check_collinearity(model.mul)
## # Check for Multicollinearity
## 
## Low Correlation
## 
##              Term  VIF Increased SE Tolerance
##         uhc_index 1.08         1.04      0.93
##  inter.travel.red 1.10         1.05      0.91
##     wear.mask.red 1.47         1.21      0.68
##   test.policy.red 1.35         1.16      0.74
##              time 1.43         1.20      0.70
# 
 pp_check(model.mul, 250)

library(glmmTMB)
## Warning in checkMatrixPackageVersion(): Package version inconsistency detected.
## TMB was built with Matrix version 1.3.4
## Current Matrix version is 1.3.3
## Please re-install 'TMB' from source using install.packages('TMB', type = 'source') or ask CRAN for a binary version of 'TMB' matching CRAN's 'Matrix' package
## 
## Attaching package: 'glmmTMB'
## The following object is masked from 'package:VGAM':
## 
##     betabinomial
model.glm <- glmmTMB(week.case.eff ~ uhc_index  + inter.travel.red  + wear.mask.red +  test.policy.red + 
                       time + (time|iso_code), data = data0.anal, family = beta_family(link="probit"), REML = TRUE)

summary(model.glm)
##  Family: beta  ( probit )
## Formula:          
## week.case.eff ~ uhc_index + inter.travel.red + wear.mask.red +  
##     test.policy.red + time + (time | iso_code)
## Data: data0.anal
## 
##      AIC      BIC   logLik deviance df.resid 
##  -2679.9  -2619.9   1354.0  -2707.9      536 
## 
## Random effects:
## 
## Conditional model:
##  Groups   Name        Variance  Std.Dev. Corr  
##  iso_code (Intercept) 0.1224981 0.3500         
##           time        0.0002592 0.0161   -0.79 
## Number of obs: 540, groups:  iso_code, 30
## 
## Dispersion parameter for beta family (): 43.3 
## 
## Conditional model:
##                    Estimate Std. Error z value Pr(>|z|)    
## (Intercept)       -0.654016   0.282376  -2.316 0.020552 *  
## uhc_index         -0.008340   0.004069  -2.050 0.040405 *  
## inter.travel.red1 -0.030012   0.061666  -0.487 0.626478    
## inter.travel.red2 -0.162898   0.067333  -2.419 0.015551 *  
## wear.mask.red1    -0.238948   0.063963  -3.736 0.000187 ***
## wear.mask.red2    -0.136099   0.076365  -1.782 0.074714 .  
## wear.mask.red3    -0.080739   0.060022  -1.345 0.178572    
## test.policy.red1  -0.060338   0.057882  -1.042 0.297213    
## test.policy.red2  -0.101022   0.068486  -1.475 0.140189    
## time              -0.034176   0.004858  -7.036 1.99e-12 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
p_load(plotly, modelbased)

data0.anal1$Predicted <- estimate_response(model.glm)$Predicted
## Warning in get_predicted.glmmTMB(model, data = data, predict = predict, : `predict = 'prediction'` is currently not available for glmmTMB models.
##   Changing to `predict = 'expectation'`.
plot2 <- data0.anal1 %>%
ggplot() +
geom_line(aes(x = log(week.case.eff), y = log(week.case.eff)), linetype = "dashed") +
geom_point(aes(x = log(week.case.eff) , y = log(Predicted), key=iso_code), color = "red") +
ylab("wADGR (predicted)") + xlab("wADGR (observed)")
## Warning: Ignoring unknown aesthetics: key
plot2

ggplotly(plot2, source = "select", tooltip = c("key") )
# Develop function marginal 
amex<- function(ame){
  x1 <- data.frame(list("AME" = ame*100, "factors" = as.character(substitute(ame))))
}
summary(model.glm)
##  Family: beta  ( probit )
## Formula:          
## week.case.eff ~ uhc_index + inter.travel.red + wear.mask.red +  
##     test.policy.red + time + (time | iso_code)
## Data: data0.anal
## 
##      AIC      BIC   logLik deviance df.resid 
##  -2679.9  -2619.9   1354.0  -2707.9      536 
## 
## Random effects:
## 
## Conditional model:
##  Groups   Name        Variance  Std.Dev. Corr  
##  iso_code (Intercept) 0.1224981 0.3500         
##           time        0.0002592 0.0161   -0.79 
## Number of obs: 540, groups:  iso_code, 30
## 
## Dispersion parameter for beta family (): 43.3 
## 
## Conditional model:
##                    Estimate Std. Error z value Pr(>|z|)    
## (Intercept)       -0.654016   0.282376  -2.316 0.020552 *  
## uhc_index         -0.008340   0.004069  -2.050 0.040405 *  
## inter.travel.red1 -0.030012   0.061666  -0.487 0.626478    
## inter.travel.red2 -0.162898   0.067333  -2.419 0.015551 *  
## wear.mask.red1    -0.238948   0.063963  -3.736 0.000187 ***
## wear.mask.red2    -0.136099   0.076365  -1.782 0.074714 .  
## wear.mask.red3    -0.080739   0.060022  -1.345 0.178572    
## test.policy.red1  -0.060338   0.057882  -1.042 0.297213    
## test.policy.red2  -0.101022   0.068486  -1.475 0.140189    
## time              -0.034176   0.004858  -7.036 1.99e-12 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
coef_model <- summary(model.glm)$coefficients$cond[,1][-1]

coef <- as.data.frame(coef_model)

# AME time
AME.time <- coef_model[9] * mean(dnorm(predict(model.glm, type="link")))
z9 <- amex(AME.time)

# uhc index
AME.uhc_idx <- coef_model[1]* mean(dnorm(predict(model.glm, type="link")))
z1 <- amex(AME.uhc_idx)

# travel
trv0 <- predict(model.glm, newdata=data.frame(list("time"=data0.anal1$time,
"uhc_index"=data0.anal1$uhc_index,
"inter.travel.red"="0",
"wear.mask.red"=data0.anal1$wear.mask.red, 
"test.policy.red"=data0.anal1$test.policy.red,
"iso_code"=data0.anal1$iso_code) ), type = "response")

trv1 <- predict(model.glm, newdata=data.frame(list("time"=data0.anal1$time,
"uhc_index"=data0.anal1$uhc_index,
"inter.travel.red"="1",
"wear.mask.red"=data0.anal1$wear.mask.red, 
"test.policy.red"=data0.anal1$test.policy.red,
"iso_code"=data0.anal1$iso_code) ), type = "response")
trv2 <- predict(model.glm, newdata=data.frame(list("time"=data0.anal1$time,
"uhc_index"=data0.anal1$uhc_index,
"inter.travel.red"="2",
"wear.mask.red"=data0.anal1$wear.mask.red, 
"test.policy.red"=data0.anal1$test.policy.red,
"iso_code"=data0.anal1$iso_code) ), type = "response")

trv10 <- trv1 - trv0 
trv20 <- trv2 - trv0 


AME.trv10<-mean(trv10)
AME.trv20<-mean(trv20)

z2<-amex(AME.trv10)
z3<-amex(AME.trv20)


#wear.mask.red 
mask0 <- predict(model.glm, newdata=data.frame(list("time"=data0.anal1$time,
"uhc_index"=data0.anal1$uhc_index,
"inter.travel.red"=data0.anal1$inter.travel.red,
"wear.mask.red"="0", 
"test.policy.red"=data0.anal1$test.policy.red,
"iso_code"=data0.anal1$iso_code) ), type = "response")
mask1 <- predict(model.glm, newdata=data.frame(list("time"=data0.anal1$time,
"uhc_index"=data0.anal1$uhc_index,
"inter.travel.red"=data0.anal1$inter.travel.red,
"wear.mask.red"="1", 
"test.policy.red"=data0.anal1$test.policy.red,
"iso_code"=data0.anal1$iso_code) ), type = "response")
mask2 <- predict(model.glm, newdata=data.frame(list("time"=data0.anal1$time,
"uhc_index"=data0.anal1$uhc_index,
"inter.travel.red"=data0.anal1$inter.travel.red,
"wear.mask.red"="2", 
"test.policy.red"=data0.anal1$test.policy.red,
"iso_code"=data0.anal1$iso_code) ), type = "response")
mask3 <- predict(model.glm, newdata=data.frame(list("time"=data0.anal1$time,
"uhc_index"=data0.anal1$uhc_index,
"inter.travel.red"=data0.anal1$inter.travel.red,
"wear.mask.red"="3", 
"test.policy.red"=data0.anal1$test.policy.red,
"iso_code"=data0.anal1$iso_code) ), type = "response")

mask10<-mask1-mask0
mask20<-mask2-mask0
mask30<-mask3-mask0

AME.mask10<-mean(mask10)
AME.mask20<-mean(mask20)
AME.mask30<-mean(mask30)


z4<-amex(AME.mask10)
z5<-amex(AME.mask20)
z6<-amex(AME.mask30)

# test 
test0<- predict(model.glm, newdata=data.frame(list("time"=data0.anal1$time,
"uhc_index"=data0.anal1$uhc_index,
"inter.travel.red"=data0.anal1$inter.travel.red,
"wear.mask.red"=data0.anal1$wear.mask.red, 
"test.policy.red"="0",
"iso_code"=data0.anal1$iso_code) ), type = "response")
test1<- predict(model.glm, newdata=data.frame(list("time"=data0.anal1$time,
"uhc_index"=data0.anal1$uhc_index,
"inter.travel.red"=data0.anal1$inter.travel.red,
"wear.mask.red"=data0.anal1$wear.mask.red, 
"test.policy.red"="1",
"iso_code"=data0.anal1$iso_code) ), type = "response")
test2<- predict(model.glm, newdata=data.frame(list("time"=data0.anal1$time,
"uhc_index"=data0.anal1$uhc_index,
"inter.travel.red"=data0.anal1$inter.travel.red,
"wear.mask.red"=data0.anal1$wear.mask.red, 
"test.policy.red"="2",
"iso_code"=data0.anal1$iso_code) ), type = "response")


test10<-test1-test0
test20<-test2-test0


AME.test10<-mean(test10)
AME.test20<-mean(test20)


z7<-amex(AME.test10)
z8<-amex(AME.test20)


z_merge <- rbind(z1,z2,z3,z4,z5,z6,z7,z8,z9)

term3 <- c(
"UHC Index\n",
"Ban arrivals from\nsome regions\n",
"Ban on all regions\n",
'Recommend or Required wear\nmask in some public spaces',
'Required wear mask in all\n public spaces with \nother people present',
'Required wear mask outside\n the home at all times',
"Testing of anyone \nshowing symptoms\n",
"Open public testing\n",
"Time\n"
  )
p_load(broom.mixed)
sum_table <- tidy(model.mul, conf.int = T)

# create new table summary of model
new_table <- sum_table[2:10,]
 
new_table1 <- cbind(new_table, term3, z_merge)

new_table2 <- new_table1 %>% select(term, term3, AME, estimate, conf.low, conf.high)

new_table2$term3 <- factor(new_table2$term3, levels = new_table2$term3)


p.vac.period<- ggplot(new_table2, aes(x = term3, y = estimate,
                     ymin = conf.low, ymax = conf.high)) +
    geom_hline( yintercept = 0, color = 'Black' ) + ylim(-0.5,0.5) + 
    geom_linerange(size =1, color = 'steelblue') + geom_point(size=1.5, color ='steelblue') + coord_flip() + 
  ylab("Probit Transformation of Average Weekly Growth Rate ") + 
  xlab("Predictors") +
  theme_minimal() + 
  theme(strip.background = element_rect(NA), 
  axis.text = element_text(size = 20),
  axis.title = element_text(size = 20),
  legend.text = element_text(size = 10),
  legend.title = element_text(size = 10),
  strip.text = element_text(size = 8),
  legend.position = "top")
p.vac.period

ame.plot <- ggplot(data = new_table2, aes(x= term3, y = AME)) + geom_bar(stat="identity",fill="steelblue") + coord_flip() +
  geom_hline( yintercept = 0, color = 'Black' ) +
  geom_text(aes(label=round(AME,2)), hjust= 1.05, color="black", size=5)+
  ylab("Average marginal effect (%)") +
  ylim(-3,3) +
  xlab("") +
  theme_minimal() + 
  theme(strip.background = element_rect(NA), 
        axis.text.y=element_blank(),
  axis.text = element_text(size = 20),
  axis.title = element_text(size = 20),
  legend.text = element_text(size = 10),
  legend.title = element_text(size = 10),
  strip.text = element_text(size = 8),
  legend.position = "top")
ame.plot

cairo_pdf(filename = "F:/STUDY HCA IN TAIWAN/0. ideas for thesis/COVID19 and NPIs/Data analysis/results/plot_sep/fig3_2020_Sep29.pdf",
          width = 20, height = 10, #inch
          onefile = TRUE, family = "Segoe UI")
merge_bar <- ggarrange(p.vac.period, ame.plot, ncol = 2, nrow = 1, widths = c(1.5,1))
merge_bar
dev.off()
## png 
##   2